%% Solve model


% Save parameters
cd   'model'
save paramfile sig_eps_i sig_eps_a betta siggma varphi eta ...
     theta_1 theta_2 theta_3 theta_4 theta_5 phi_pi phi_y pibar ...
     rho_a rho_i pbar_1 pbar_2 pbar_3 pbar_4 pbar_5 ...
     ybar mctildebar pstarbar_1 pstarbar_2 pstarbar_3 pstarbar_4 pstarbar_5 ...
     sbar_1 sbar_2 sbar_3 sbar_4 sbar_5 sbar ...
     markupbar_1 markupbar_2 markupbar_3 markupbar_4 markupbar_5 ...
     psibar_1 psibar_2 psibar_3 psibar_4 psibar_5 ...
     phibar_1 phibar_2 phibar_3 phibar_4 phibar_5 ...
     abar nbar tfpbar ibar i_annbar pi_annbar posmp


% Solve 3rd order approximation around deterministic steady state
order = 3;
dynare model_03_rotemberg.mod nograph -Dalternative_taylor_rule=0
cd ..


% Compute linearized model dynamics around stochastic steady state
type = 'stochastic_steady_state';
[deflect_sss] = compute_deflected_linear_approximation(M_,options_,oo_,type);
[deflect_sss] = compute_deflected_irf(M_,options_,deflect_sss,var_list_,0); % last argument: do plots


% Position of variables
fun_var_pos


% IRFs: MP shock
irf_mp_sss  = deflect_sss.irf(1:M_.endo_nbr,:)';

% IRFs: tech shock
irf_tfp_sss = deflect_sss.irf(M_.endo_nbr+1:M_.endo_nbr*2,:)';
fac_tfp     = sqrt(M_.Sigma_e(2,2)); 
irf_tfp_sss = irf_tfp_sss*fac_tfp;   

% Markup dispersion IRF: MP shock
pos.markupdisp    = M_.endo_nbr+1;
markup1_log_irf   = deflect_sss.y(pos.markup_1) + irf_mp_sss(:,pos.markup_1); 
markup2_log_irf   = deflect_sss.y(pos.markup_2) + irf_mp_sss(:,pos.markup_2); 
markup3_log_irf   = deflect_sss.y(pos.markup_3) + irf_mp_sss(:,pos.markup_3); 
markup4_log_irf   = deflect_sss.y(pos.markup_4) + irf_mp_sss(:,pos.markup_4); 
markup5_log_irf   = deflect_sss.y(pos.markup_5) + irf_mp_sss(:,pos.markup_5); 
markupAvg_log_irf = 1/5*markup1_log_irf + 1/5*markup2_log_irf + 1/5*markup3_log_irf + 1/5*markup4_log_irf + 1/5*markup5_log_irf;
markupdisp_irf    = 1/5*(markup1_log_irf - markupAvg_log_irf).^2 + 1/5*(markup2_log_irf - markupAvg_log_irf).^2 + 1/5*(markup3_log_irf - markupAvg_log_irf).^2 + 1/5*(markup4_log_irf - markupAvg_log_irf).^2 + 1/5*(markup5_log_irf - markupAvg_log_irf).^2;
markupAvg_sss     = 1/5*deflect_sss.y(pos.markup_1) + 1/5*deflect_sss.y(pos.markup_2) + 1/5*deflect_sss.y(pos.markup_3) + 1/5*deflect_sss.y(pos.markup_4) + 1/5*deflect_sss.y(pos.markup_5);
markupdisp_sss    = 1/5*(deflect_sss.y(pos.markup_1) - markupAvg_sss)^2 + 1/5*(deflect_sss.y(pos.markup_2) - markupAvg_sss)^2 + 1/5*(deflect_sss.y(pos.markup_3) - markupAvg_sss)^2 + 1/5*(deflect_sss.y(pos.markup_4) - markupAvg_sss)^2 + 1/5*(deflect_sss.y(pos.markup_5) - markupAvg_sss)^2;
irf_mp_sss(:,pos.markupdisp) = markupdisp_irf - markupdisp_sss;

% Markup dispersion IRF: tech shock
markup1_log_irf   = deflect_sss.y(pos.markup_1) + irf_tfp_sss(:,pos.markup_1); 
markup2_log_irf   = deflect_sss.y(pos.markup_2) + irf_tfp_sss(:,pos.markup_2); 
markup3_log_irf   = deflect_sss.y(pos.markup_3) + irf_tfp_sss(:,pos.markup_3); 
markup4_log_irf   = deflect_sss.y(pos.markup_4) + irf_tfp_sss(:,pos.markup_4); 
markup5_log_irf   = deflect_sss.y(pos.markup_5) + irf_tfp_sss(:,pos.markup_5); 
markupAvg_log_irf = 1/5*markup1_log_irf + 1/5*markup2_log_irf + 1/5*markup3_log_irf + 1/5*markup4_log_irf + 1/5*markup5_log_irf;
markupdisp_irf    = 1/5*(markup1_log_irf - markupAvg_log_irf).^2 + 1/5*(markup2_log_irf - markupAvg_log_irf).^2 + 1/5*(markup3_log_irf - markupAvg_log_irf).^2 + 1/5*(markup4_log_irf - markupAvg_log_irf).^2 + 1/5*(markup5_log_irf - markupAvg_log_irf).^2;
markupAvg_sss     = 1/5*deflect_sss.y(pos.markup_1) + 1/5*deflect_sss.y(pos.markup_2) + 1/5*deflect_sss.y(pos.markup_3) + 1/5*deflect_sss.y(pos.markup_4) + 1/5*deflect_sss.y(pos.markup_5);
markupdisp_sss    = 1/5*(deflect_sss.y(pos.markup_1) - markupAvg_sss)^2 + 1/5*(deflect_sss.y(pos.markup_2) - markupAvg_sss)^2 + 1/5*(deflect_sss.y(pos.markup_3) - markupAvg_sss)^2 + 1/5*(deflect_sss.y(pos.markup_4) - markupAvg_sss)^2 + 1/5*(deflect_sss.y(pos.markup_5) - markupAvg_sss)^2;
irf_tfp_sss(:,pos.markupdisp) = markupdisp_irf - markupdisp_sss;

% Aggregate markup IRF: MP shock
pos.agg_markup    = M_.endo_nbr+2;
irf_mp_sss(:,pos.agg_markup)  = -irf_mp_sss(:,pos.mctilde)+irf_mp_sss(:,pos.tfp);

% Aggregate markup IRF: tech shock
irf_tfp_sss(:,pos.agg_markup) = -irf_tfp_sss(:,pos.mctilde)+irf_tfp_sss(:,pos.tfp);


